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Abstract. A new anisotropic mesh adaptation strategy for finite element solution 
of elliptic differential equations is presented. It generates anisotropic adaptive meshes 
as quasi-uniform ones in some metric space, with the metric tensor being computed 
based on a posteriori error estimates proposed in j36| . The new metric tensor explores 
more comprehensive information of anisotropy for the true solution than those exist- 
ing ones. Numerical results show that this approach can be successfully applied to 
deal with poisson and steady convection-dominated problems. The superior accuracy 
and efficiency of the new metric tensor to others is illustrated on various numerical 
examples of complex two-dimensional simulations. 
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1 Introduction 

Nowadays, many computational simulations of partial differential equations (PDEs) in- 
volve adaptive triangulations. Mesh adaptation aims at improving the efficiency and the 
accuracy of numerical solutions by concentrating more nodes in regions of large solution 
variations than in other regions of the computational domain. As a consequence, the 
number of mesh nodes required to achieve a given accuracy can be dramatically reduced 
thus resulting in a reduction of the computational cost. Traditionally, isotropic mesh 
adaptation has received much attention, where almost regular mesh elements are only 
adjusted in size based on an error estimate. However, in regions of large solution gradient, 
adaptive isotropic meshes usually contain too many elements. Moreover, if the problem 
at hand exhibits strong anisotropic features that their solutions change more significantly 
in one direction than the others, like boundary layers, shock waves, interfaces, and edge 
singularities, etc.. In such cases it is advantageous to refiect this anisotropy in the dis- 
cretization by using meshes with anisotropic elements (sometimes also called elongated 



*LSEC, ICMSEC, Academy of Mathematics and Systems Science, CAS, Beijing 100080, China email: 
hhxie@lsec.cc.ac.cn 

^Department of Mathematics, Central China Normal University, Wuhan 430079, China email: 
yinxb@lsec.cc.ac.cn 



1 



elements). These elements have a small mesh size in the direction of the rapid variation 
of the solution and a larger mesh size in the perpendicular direction. Indeed anisotropic 
meshes have been used successfully in many areas, for example in singular perturbation 
and flow problems [IlEllllEOlEIlEIlEH] and in adaptive procedures I5l|6l|8l[l0l[22l|3ll|32]. 
For problems with very different length scales in different spatial directions, long and thin 
triangles turn out to be better choices than shape regular ones if they are properly used. 

Compared to traditionally used isotropic ones, anisotropic meshes are more difficult to 
generate, requiring a full control of both the shape, size, and orientation of elements. It 
is necessary to have as much information as possible on the nature and local behavior of 
the solution. We need to convert somehow what we know about the solution to something 
having the dimension of a length and containing directional information. In practice, 
they are commonly generated as quasi-uniform meshes in the metric space determined 
by a tensor (or a matrix-valued function) specifying the shape, size, and orientation of 
elements on the entire physical domain. 

Such a metric tensor is often given on a background mesh, either prescribed by the user 
or chosen as the mesh from the previous iteration in an adaptive solver. So far, several 
meshing strategies have been developed for generating anisotropic meshes according to a 
metric prescription. Examples include blue refinement |27ll28j . directional refinement [32], 
Delaunay-type triangulation method [5l[6l[10l[31], advancing front method pS|, bubble 
packing method [33], local refinement and modification [20^33]. On the other hand, vari- 
ational methods have received much attention in the recent years typically for generating 
structured meshes as they are especially well suited for finite difference schemes. In these 
approaches, an adaptive mesh is considered as the image of a computational mesh under 
a coordinate transformation determined by a functional [71[T51[231[2S1[2S1[2S] . Readers are 
referred to [17] for an overview. 

Among these meshing strategies, the definition of the metric tensor based on the Hessian 
of the solution seems nowadays commonly generalized in the meshing community. This 
choice is largely motivated by the interesting numerical results obtained by the results of 
D'Azevedo [13], D'Azevedo and Simpson |14) on linear interpolation for quadratic functions 
on triangles. For example, Castro-Dzaz et al. [TO], Habashi et al. [20], and Remade et 
al. [23] define their metric tensor as 



To guarantee its positive definiteness and avoid unrealistic metric, A4 in (jl.ip is modified 
by imposing the maximal and minimal edge lengths. Hecht [22] uses 




(1.1) 



where the Hessian of function u has the eigen-decomposition H{u) = R diag(Ai, A2)-R^. 
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for the relative error and 

M= ' , l^"'!, , (1.3) 

eo • Coef sup(n) - mi(u) 

for the absolute error, where eo, Coef, and CutOff are the user specified parameters used 
for setting the level of the linear interpolation error (with default value 10~^), the value 
of a multiplicative coefficient on the mesh size (with default value 1), and the limit value 
of the relative error evaluation (with default value 10^^), respectively. In [I9|, George and 
Hecht define the metric tensor for various norms of the interpolation error as 

M={^Y4 (1.4) 



eo/ \ \\\ ^ 

where cq is a constant, eo is a given error threshold, and v = 1 for the norm and the 
semi-norm and v = 1/2 for norm of the error. It is emphasized that the definitions 
()l.ip - (|1.3p are based on the results of [13] while ()1.4p largely on heuristic considerations. 
Huang [24] develop the metric tensors as 



i.(«)det(/ + l 

for the LP' norm and 



1 



A4o = 3- (-)det(/ + ^|i/(u)|) ' I + ^\H{u)\ (1.5) 
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Mi = --[-] I+-\H{u)\ (1.6) 



for the semi-norm. 

This list is certainly incomplete, but from the papers we can draw two conclusions. 
First, compared with isotropic mesh, significant improvements in accuracy and efficiency 
can be gained when a properly chosen anisotropic mesh is used in the numerical solution 
for a large class of problems which exhibit anisotropic solution features. Second, there are 
still challenges to access fully anisotropic information from the computed solution. 

The objective of this paper is to develop a new way to get metric tensors for anisotropic 
mesh generation in two dimension, which explores more comprehensive anisotropic in- 
formation than some exist methods. The development is based on the error estimates 
obtained in our recent work [36j on linear interpolation for quadratic functions on trian- 
gles. These estimates are anisotropic in the sense that they allow a full control of the 
shape of elements when used within a mesh generation strategy. 

The application of the error estimates of |36] to formulate the metric tensor, is based 
on two factors: on the one hand, as a common practice in the existing anisotropic mesh 
generation codes, we assume that the anisotropic mesh is generated as a quasi-uniform 
mesh in the metric tensor A^, i.e., a mesh where the elements are equilateral or isosceles 
right triangle or other quasi-uniform triangles (shape requirement) in M and unitary in 
size (size requirement) in AA. On the other hand, the anisotropic mesh is required to 
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minimize the error for a given number of mesh points (or equidistribute the error). Then 
is constructed from these requirements. We will establish new metric tensors as 



N 



■ det H 



n 



for the norm and 
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for the semi-norm, where N is the number of elements in the triangulation. Under 
the condition Ti = I + ^\H{u)\, our metric tensor A^o for the norm is similar to (jl.Sp . 
However, the new metric tensor TWi is essentially different with those metric tensors 
mentioned above. The difference lies on the term 



Vdet(^) 



(1.9) 



which indicates our metric tensor explores more comprehensive anisotropic information of 
the solution when the term p.9p varies significantly at different points or elements. In 
addition, numerical results will show that the more anisotropic the solution is, the more 
obvious the superiority of the new metric tensor is. 

The paper is organized as follows. In Section 2, we describe the anisotropic error 
estimates on linear interpolation for quadratic functions on triangles obtained in the recent 
work [36j. The formulation of the metric tensor is developed in Section 3. Numerical results 
are presented in Section 4 for various examples of complex two-dimensional simulations. 
Finally, conclusions are drawn in Section 5. 



2 Anisotropic estimates for interpolation error 

Needless to say, the interpolation error depends on the solution and the size and shape 
of the elements in the mesh. Understanding this relation is crucial for the generation 
of efficient and effective meshes for the finite element method. In the mesh generation 
community, this relation is studied more closely for the model problem of interpolating 
quadratic functions. This model is a reasonable simplification of the cases involving general 
functions, since quadratic functions are the leading terms in the local expansion of the 
linear interpolation errors. For instance, Nadler |30) derived an exact expression for the 
L^-norm of the linear interpolation error in terms of the three sides £i, £2, and £3 of the 
triangle K (see Figure [1]) , 



\U - Ul\\r2 



2 ^2 



hiK) = Y^[[di + d2 + d3) +did2 + d2d3 + did3\. (2.1) 
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where \K\ is the area of the triangle, di = £i ■ Hii with H being the Hessian of u. Bank 
and Smith [3] gave a formula for the iif^-seminorm of the linear interpolation error 

||V(n-n7)||^2(^) = ^v-5v, 

where v = [di, ^2, c^s]^, 

2^1-^3 2^2-^3 + + 

Cao derived two exact formulas for the //^-seminorm and L^-norm of the interpolation 
error in terms of the area, aspect ratio, orientation, and internal angles of the triangle. 

Chen, Sun and Xu [11] obtained the error estimate 

||V(u-n/)||?p.o^ < CiV" nil {/detail pn_ ,l<j9<oo, 

where the constant C does not depend on u and N. They also showed the estimate is 
optimal in the sense that it is a lower bound if u is strictly convex or concave. 

Assuming u = Xix^ + A2y^, D'Azevedo and Simpson |12| derived the exact formula for 
the maximum norm of the interpolation error 

^||2 D12D23D31 
"^^-^^)"^-W=16AiA2|i^r 

where Dij = ii •diag(Ai, A2)^i. Based on the geometric interpretation of this formula, they 
proved that for a fixed area the optimal triangle, which produces the smallest maximum 
interpolation error, is the one obtained by compressing an equilateral triangle by factors 
\/A]" and \/Ai along the two eigenvectors of the Hessian of u. Furthermore, the optimal 
incidence for a given set of interpolation points is the Delaunay triangulation based on 
the stretching map (by factors ^/Xi and ^/Xi along the two eigenvector directions) of the 
grid points. Rippa [S] showed that the mesh obtained this way is also optimal for the 
L^-norm of the error for any 1 < p < 00. 

The element-wise error estimates in the following theorem are developed in using 
the theory of interpolation for finite elements and proper numerical quadrature formula. 

Theorem 2.1. Letu be a quadratic function anduj is the Lagrangian linear finite element 
interpolation of u. The following relationship holds: 



|V(n - uj)\\l,^^^ = ^(£,+1 • Hi,+2W, 



' ' i=l 

where we prescribe i + 3 = i,i — 3 = i. 



5 



Naturally, 



VI 



^ 48|K| ^ 



• HI 



i+2j 



is set as the a posteriori estimator for ||V(u — uj)\\i^2i^Qy Numerical experiments in [36] 
show that the estimators are always asymptotically exact under various isotropic and 
anisotropic meshes. 



3 Metric tensors for anisotropic mesh adaptation 

We now use the results of Theorem 12.11 to develop a formula for the metric tensor. As 
a common practice in anisotropic mesh generation, we assume that the metric tensor, 
M (x) , is used in a meshing strategy in such a way that an anisotropic mesh is generated 
as a quasi-uniform mesh in the metric determined by A^(x). Mathematically, this can be 
interpreted as the shape, size and equidistribution requirements as follows. 

The shape requirement. The elements of the new mesh, Th, are (or are close to 
being) equilateral in the metric. 

The size requirement. The elements of the new mesh Th have a unitary volume in 
the metric, i.e., 

/ y'det{M{x))dx =1, WK G Th- (3.1) 
Jk 

The equidistribution requirement. The anisotropic mesh is required to minimize 
the error for a given number of mesh points (or equidistribute the error on every element). 
We now state the most significant contributions of this paper. 

3.1 Metric tensor for the semi-norm 

Assume that H{x) is a symmetric positive definite matrix on every point x. Let A^i(x) = 
6iMi(x) where Oi is to be determined. Here Mi(x) is often called the monitor function. 
Both the monitor function Mi(x) and metric tensor A^i(x) play the same role in mesh 
generation, i.e., they are used to specify the size, shape, and orientation of mesh elements 
throughout the physical domain. The only difference lies in the way they specify the size of 
elements. Indeed, Mi(x) specifies the element size through the equidistribution condition, 
while 7Wi(x) determines the element size through the unitary volume requirement ()3.ip . 

Assume that H{x) is a constant matrix on K, denoted by Hk, then so does Mi(x), 
denoted by Mi^k- Since Hk is a symmetric positive definite matrix, we consider the 
singular value decomposition Hk = R^AR, where A = diag(Ai, A2) is the diagonal matrix 
of the corresponding eigenvalues (Ai,A2 > 0) and R is the orthogonal matrix having as 
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Figure 1: Affine map x = J^^-x from K to the reference triangle K. 



rows the eigenvectors of Hk- Denote by Fk and tx the matrix and the vector defining the 
invertible affine map x = J^/^(x) = -Fftrx + from the generic element K to the reference 
triangle K (see Figured]). 

1 

Let Mi^K = CkHk, then Fx = C^-Aai? and Mij^ = FJ^Fk- Mathematically, the 
shape requirement can be expressed as 

= L and cos^j = ^^+'^^^^+'^ ^ ^ ^ 1,2,3, 
where L is a constant for every element K. Together with Theorem 12.11 we have 

3 



r4 3 



mK\c^ .^^ 



L4 



32^/3C|, 

To satisfy the equidistribution requirement, let 



+ \/A2/Ai). 



|V(n-^z/)||i.(^) = ( J] e|,)/7V = 6?/iV, 
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where is the number of elements of Th- Then 

r2/ N 



Ck = 



32V3ef 



C 



tr(H] 



Ly/d^tCff) 



where C is a constant which depends on the prescribed error and the numbers of elements. 
Generally speaking, -ff(x) can not be a constant matrix on K , Mi(x) should be the form 



Mi(x) 



r tiiH) 1 


1 

2 


Hiu) 




Lvdet(/7)J 









VAi(x)/A2(x) + v/A2(x)/Ai(x 



H(u) 



since Mi(x) can be modified by multiplying a constant. Note that Ai(x) and A2(x) are 
the corresponding eigenvalue of H{u) at point x. 

To establish 7Wi(x), the size requirement (|3.ip should be used, which leads to 

Oi / /9i(x)fix = 1, 
Jk 

where 

pi(x) = Vdet(Mi(x)). 
Summing the above equation over all the elements of Th, one gets 

9iai = N, 



where 



Thus, we get 



and as a consequence, 



cji = / pi(x)(ix. 
Jn 



N 
0-1 



A1i(x) 



N 



tr(H] 



H{u) 



N 



VAi(x)/A2(x) + VA2(x)/Ai( 



(y\ LydetpO 
3.2 Metric tensor for the L? norm 



Using the expression ()2.ip for the L^-norm of the linear interpolation error derived by 
Nadler [30] and the condition (13.111. we have 



W-'^ifb{K) = "^[(X^^*) + did2 + d2d3 + dids 

i=l 



180 
\K 



180Ci 



15C|. ~ ISCfTCT ~ 60C|VAIA^' 
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To satisfy the equidistribution requirement, let 

\W-Ml^K) = { E el)/N = el/N. 
Ken 

Using similar argument with last subsection, we easily get the monitor function 



and the metric tensor 



for the norm. 



Mo(x) = det{H 



N 



Hiu] 



Xo(x) = — det i7 H{u) 



3.3 Practice use of metric tensors 

So far we assume that i^(x) is a symmetric positive definite matrix at every point. However 
this assumption doesn't hold in many cases. In order to obtain a symmetric positive 
definite matrix, the following procedure are often implemented. First, the Hessian H is 
modified into \H\ = diag( | Ai |, IA2I )ii by taking the absolute value of its eigenvalues 
( |2I])- Since \H\ is only semi-positive definite, A^o and M.i cannot be directly applied to 
generate the anisotropic meshes. To avoid this difficulty, we regularize the expression with 
two flooring parameters ao > for A^o and ai > for A^i, respectively ( [23]). Replace 
\H\ with 

n = ail + \H\, 
then we get the modified metric tensor 



A^i(x) 



(Ti L^det(^) 



n 



(3.2) 



Similarly, replacing \H\ with 



leads to 



n = aoI+\H\ 



N 

7Wo(x) = —det('H 
(^0 



n 



(3.3) 



The two modified metric tensors ()3.2p and (j3.3p are suitable for practical mesh generation. 
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4 Numerical experiments 



We have elaborated that our metric tensor explores more anisotropic information for a 
given problem. Then it is crucial to check that whether our metric tensor exhibits better 
than those without the term 



tr(n) 
7dit(?7) 



or VAi(x)/A2(x) + VA2(x)/Ai(x) 



(4.1) 



4.1 Mesh adaptation tool 

All the presented experiments are performed by using the BAMG software [22]. Given a 
background mesh and an approximation solution, BAMG generates the mesh according 
to the metric. The code allows the user to supply his/her own metric tensor defined on a 
background mesh. In our computation, the background mesh has been taken as the most 
recent mesh available. 



4.2 Comparisons between two types of metric tensors 

Specifically, in every example, the PDE is discretized using linear triangle finite elements. 
Two serials of adaptive meshes of almost the same number of elements are shown (the 
iterative procedure for solving PDEs is shown in Figure [2]) . The finite element solutions 











Compute metric 




Generate new mesh 


Given a mesh 


> 


Solve PDE 




tensor on the 
current mesh 




according to the 
metric tensor 




> 


> 



iteration 



Figure 2: Iterative procedure for adaptive mesh solution of PDEs. 



are computed by using the metric tensor of modified Hessian 



N 
0-1 



\H\ 



and the new metric tensor 



-M?(x) 



0"! 



Vdetpl) 



H\ 



(4.2) 



(4.3) 



Notice that the formulas of the metric tensors involve second order derivatives of the 
physical solution. Generally speaking, one can assume that the nodal values of the solution 
or their approximations are available, as typically in the numerical solution of partial 
differential equations. Then, one can get approximations for the second order derivatives 
using gradient recovery techniques (such as 00] and [35]) twice or a Hessian recovery 
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technique (using piecewise quadratic polynomials fitting in least-squares sense to nodal 
values of the computed solution [37]) just once, although their convergence has been 
analyzed only on isotropic meshes. In our computations, we use the technique [IQ] twice. 

In the current computation, each run is stopped after 10 iterations to guarantee that 
the adaptive procedure tends towards stability. 

Denote by nbt the number of the elements (triangles in 2D case) in the current mesh. 
The number of triangles or nodes is adjusted when necessary by trial and errors through 
the modification of the multiplicative coefficient of the metric tensors. 

Example 4.1. This example, though numerically solved in 17 = (0, 1) x (0, 1), is in fact 
one-dimensional: 

-kAu + -— = f, 
ox I 

with / = 0, u{xi = 0,^2) = 0, u{xi = 1,^2) = 1, and |^ = along the top and bottom 
sides of the square (taken from [8]). The exact solution is given by 

1 - e « 



1 ! 

1 - e« 

with a boundary layer of width k at xi ~ 1. We have set n = 0.0015, so that convection is 
dominant and the Galerkin method yields oscillatory solutions unless the mesh is highly 
refined at the boundary layer. 

Figure [3] (a) compares semi-norms of the discretization errors for the two metric 
tensors. In (b) semi-norm of error is computed by the difference between the Hessian 
of u and the recovered one, which exhibits the quality of meshes to some certain extent. 
Example 4.2. We study the Poisson equation (taken from jl6] ) 
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Figure 4: Example 4.2. (a) Anisotropic mesh obtained with the metric tensor A4™(x): 
nbt = 4244, |e|j:^i = 0.3727, and \e\jj2 = 1762. (b) Anisotropic mesh obtained the metric 
tensor M^i^): nbt = 4243, |e|^^i = 0.2842, and \e\H2 = 1101. 

( -Au = f, in J7 EE (0, 1) X (0,1), 
] u = 0, on dQ, 

where / has been chosen in such a way that the exact solution u(x) = [1 — e~"^^ — (1 — 
e~")xi]4x2(l — X2) and a is chosen to be 1000. Notice that the function u exhibits an 
exponential layer along the boundary xi = with an initial step of a. Figure H] contains 

(a) (b) 



10* 




10^ 10' 10' 10^ 10^ 10^ 10* 10^ 

nbt (Numbers of Elements) nbt (Numbers of Elements) 



Figure 5: Example 4.2. The and semi-norms of discretization error are plotted as 
functions of the number of elements {nbt) in (a) and (b), respectively. 

two meshes obtained by the two different metric tensors (j4.2p and (j4.3p . It is easily 
seen that the two meshes are obvious different. The former explores more anisotropic 
features of the solution u than the latter. In other words, the term ()4.ip complements 
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Figure 6: Example 4.3. /3 = 40, (a) Anisotropic mesh obtained with the metric tensor 
Mfix): nbt = 892, |e|j^i = 0.2581, |e|^^2 = 102.0. (b) Anisotropic mesh obtained with 
the metric tensor A^^(x): nbt = 891, |e|//i = 0.1893, |e|jy2 = 57.57. 

more comprehensive information of the exact solution. Figure [S] contains and 
semi-norms of error similar to Figure [3l 

Example 4.3. Let Q = (0, 1) x (0, 1), and u be the solution of 

-An = (3{l3 - l)xf~^(l - xf) + 2/3(2/3 - l)xf~^{l - ) 

with boundary conditions u = along the sides xi = 1 and X2 = 1, and du/dn = 
along xi = and X2 = (taken from [8]). The exact solution n(x) = (1 — x^)(l — xlf) 
exhibits two boundary layers along the right and top sides, the latter being stronger than 
the former. 

Figure E] contains two meshes obtained by the two different metric tensors from solution 
perspective, where the difference is obvious. This comparison, again, indicates our metric 
tensor explores more comprehensive anisotropic information of the solution u. Figure[7]and 
Figure [8] contain and semi-norms of error, respectively, under various parameters 
/3. We know that the larger the /3 becomes, the more significant the anisotropy of the 
solution exihibits. Then we can conclude that the more anisotropic the solution is, the 
more obvious the difference is. In fact, the difference lies on the term (|4.ip which indicates 
our metric tensor explores more comprehensive anisotropic information of the solution 
when the term varies significantly at different points or elements. 



5 Conclusions 



In the previous sections we have developed a new anisotropic mesh adaptation strategy for 
finite element solution of elliptic differential equations. Note that the new metric tensor 
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(a) (b) 




10^ 10' 10* 10 lo'' 10° 10' 10^ 

nbt (Numbers of Elements) nbt (Numbers of Elements) 



Figure 7: Example 4.3. The H semi-norms of discretization error are plotted as functions 
of nbt under the conditions (a) /3 = 5, (b) /3 = 10, (c) f3 = 20, (d) f3 = 40, respectively. 

A^o for the norm is similar to (jl.Sp ( |24j). However, the new metric tensor A^i is 
essentially different with those metric tensors existed. The difference lies on the term (jl.Op 
which indicates our metric tensor explores more comprehensive anisotropic information of 
the solution when the term varies significantly at different points or elements. Numerical 
results also show that this approach is superior to the existing ones to deal with poisson 
and steady convection-dominated problems. 
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